#!/usr/bin/python

# python version of semi-empirical element model ASMITA 
# (Aggregated Scale Morphological Interaction between a 
# Tidal inlet and the Adjacent coast, Stive et al. 1998)

# python version of asmita code by Hannes Rennau
# rewritten from fortran code by Wang et al. and
# modified input file handling, ...


import time
import numpy as numpy
import os,sys
import netCDF4
from scipy.interpolate import interp1d
 
def asmita(dt,h,t0,t1,a,slr,cce,rn,p,tyear):
    # hannes: b is not known in function asmita of fortran code
    #         and rhs not known in lubksb
    # rhs wird in lubksb gelesen, aber als b verarbeitet, nach
    # funktionsaufruf is dann aber trotzdem rhs modifiziert ....
    global h0,dlta,delta,nmax,amat,indx,ve,v,vv,c0,tran,b,rhs,c,trsp,ce,w,q,en,ceglob,vcum,s
    global ifrom,jto,alfa,calls,trange

    # exchange over three elements
    ijmax=3

    # integer
    nt=1
    n=nmax
    for i in range(1,n+1):
        c0[i]=0.0
        ceglob[i]=cce
        ce[i]=cce
        en[i]=rn
        s[i]=(-1)**(i+1)
        c[i]=ceglob[i]
        w[i]=ws[i]
        # all vars same as fortran code (print c0[i],ceglob[i], ce[i],en[i], s[i],c[i], w[i])


    for i in range(0,n+1):
        for j in range(0,n+1):
            delta[i,j]=0.0
            q[i,j]=0.0
            # all vars same as in fortran code (print 'OUT: ',delta[i,j],q[i,j])
    
    nc0=1
    nci=2
    nctot=nc0+nci
    for i in range(1,nc0+1):
        j=3
        c0[j]=cce
        ifrom[i]=0
        jto[i]=j

    for i in range(1,nc0+1):
        j=3
        delta[0,j]=dlta[3]
        q[0,j]=0.

    i=1
    j=2
    delta[i,j]=dlta[1]
    q[i,j]=0.
    ifrom[nc0+1]=i
    jto[nc0+1]=j
    i=2
    j=3
    delta[i,j]=dlta[2]
    q[i,j]=0.
    ifrom[nc0+2]=i
    jto[nc0+2]=j

    for i in range(0,n): # equals fortran: do i=0,n-1
        j=i+1
        while j<n+1:
            q[i,j]=-q[i,j]*86400.
            q[j,i]=-q[i,j]
            # this is effectively sign(1.,i-j)*Q(i,j)
            # there are 86400 seconds in a day
            delta[i,j]=delta[i,j]*86400.
            delta[j,i]=delta[i,j]
            j+=1

    # dt = time step in days
    # nt = number of time steps
    dt=(t1-t0)*360.
    if dt<=3.:
        nt=1
    else:
        nt=int(dt/3.)+1
        dt=dt/nt

    # fill the matrix to solve c from ce (amat*c = rhs)
    for i in range(1,n+1):
        for j in range(1,n+1):
            amat[i,j]=q[i,j]+delta[i,j]
        amat[i,i]=-delta[0,i]-w[i]*a[i]
        for k in range(1,n+1):
            amat[i,i]=amat[i,i]-delta[i,k]

    # LU decomposition of the matrix
    ludcmp(n,nmax)

    # start time loop #
    for jt in range(1,nt+1):
        # calculate equilibrium volumes
        mkve(p)

        # zeeget (zeespiegelrijzing):
        dh=(trange-h0)/nt
        v[1]=v[1] - a[1] * (slr/nt) + a[1]*dh/2
#        print '1',v[2]
        v[2]=v[2] + a[2] * (slr/nt)
        v[3]=v[3] - a[3] * (slr/nt)

        # equilibrium concentrations
        for j in range(1,n+1):
            ce[j]=ceglob[j]*(v[j]/ve[j])**(s[j]*en[j]) #alculate local equilibroum conc. ce, eq. (14) 

        # right hand side of conservation equation
        for i in range(1,n+1):
            rhs[i]=-q[i,0]*c0[i]-delta[0,i]*c0[i]-w[i]*a[i]*ce[i]

        # concentrations
        lubksb(n,nmax)
        cpy(n)

        # forward Euler integration.
        for k in range(1,n+1):
            v[k]=v[k]+w[k]*a[k]*s[k]*(c[k]-ce[k])*dt
 
        # processing of the landsubsidance to the model results
        if lsub_yes: botd(tyear,nmax,n,t0,dt,s)
 
        transp(c0,c,q,delta,nmax,ijmax,nc0,nci,ifrom,jto)

        trsp[1]=tran[2]*1e-6
        trsp[2]=tran[3]*1e-6
        trsp[3]=-tran[1]*1e-6
#    print 'check_values: ',v[1],v[2],v[3],a[1],a[2],a[3],slr,nt,c0,ce,en,w
    # end time loop

def mkve(vtidal):
    # expressions for equilibrium volumes
    # integer nmax
    #real h,a(nmax),v(nmax),ve(nmax), alfa(nmax)
    #real alfa1,alfa2,alfa3,abasin,vtidal
    global ve,alfa,a,v
    abasin=a[1]+a[2]
    vtidal=abasin*trange-v[1]

    ve[1]=alfa[1]*abasin*trange
    ve[2]=alfa[2]*vtidal**1.55 
    ve[3]=alfa[3]*vtidal**1.23

def cpy(n):
    # integer nmax,j
    # real a(nmax),b(nmax)
    #a=numpy.zeros((nmax)) # comment: definitions really necessary?
    #b=numpy.zeros((nmax))
    global c,rhs 
    for j in range(1,n+1):
        c[j]=rhs[j]

def ludcmp(n,np):
    #************************************************************************
    #* LU decomposition using partial pivoting (row permutations) is used.  *
    #*                                                                      *
    #* Given an nxn matrix a, with physical dimension np, this subroutine   *
    #* replaces it bu the LU decomposition or the rowwise permutation of    *
    #* itself. a and n are input. a is output. indx is an output vector     *
    #* which records the row permutation effected by the partial pivoting.  *
    #* d is output as +/- 1 depending on whether the number of row          *
    #* interchanges was even or odd, respectively. This routine is used in  *
    #* combination with LUBKSB to solve linear equations or invert a matrix.*
    #*                                                                      *
    #******************** numerical recipes par 2.3 *************************

    #INTEGER n,np,indx(n),NMAX
    #REAL d,a(np,np),TINY
    #PARAMETER (NMAX=500,TINY=1.0e-20)
    #INTEGER i,imax,j,k
    #REAL aamax,dum,sum,vv(NMAX)
    global nmax,amat,indx,vv
    # vv=numpy.zeros((nmax+1),dtype=('f8')) #  
    d=1
    TINY=1.0e-20
    NMAX=500
    for i in range(1,n+1):
        aamax=0.
        for j in range(1,n+1):
            if abs(amat[i,j])>aamax: aamax=abs(amat[i,j])
        #hp if (aamax.eq.0.) pause 'singular matrix in ludcmp'
        if (abs(aamax)<1.0e-30):
            time.sleep(1)
            print 'singular matrix in ludcmp'
        vv[i]=1./aamax

    for j in range(1,n+1):
        for i in range(1,j):
            sum=amat[i,j]
            for k in range(1,i):
                sum=sum-amat[i,k]*amat[k,j]
            amat[i,j]=sum
        aamax=0.
        for i in range(j,n+1):
            sum=amat[i,j]
            for k in range(1,j):
                sum=sum-amat[i,k]*amat[k,j]

            amat[i,j]=sum
            dum=vv[i]*abs(sum)
            if dum>=aamax:
                imax=i
                aamax=dum

        if j!=imax:
            for k in range(1,n+1):
                dum=amat[imax,k]
                amat[imax,k]=amat[j,k]
                amat[j,k]=dum
            d=-d
            vv[imax]=vv[j]

        indx[j]=imax

        #hp:  if(a(j,j).eq.0.)a(j,j)=TINY
        if abs(amat[j,j])<1.0e-30: amat[j,j]=TINY
        if j!=n:
            dum=1./amat[j,j]
            for i in range(j+1,n+1):
                amat[i,j]=amat[i,j]*dum

    # amat same as in fortran code hannes
    #  (C) Copr. 1986-92 Numerical Recipes Software Ea%].3.


def lubksb(n,np):
    #********************************************************************
    #*   Back substitution                                              *
    #*                                                                      *
    #* Solves the set of equations a x = b. Here a is input, not as a       *
    #* matrix but rather as its LU decomposition, determined by the         *
    #* subroutine LUDCMP. indx is input as the permutation vector returned  *
    #* by LUDCMP. b is input as the right hand-side vector b, and returns   *
    #* with the solution vector x. a, n, np and indx are not modified by    *
    #* this subroutine and can be left in place for succesive calls with    *
    #* different right-hand sides b. This routine takes into account the    *
    #* possibility that b will begin with many zero elements, so it is      *
    #* efficient for use in matrix inversion.                               *
    #******************** numerical recipes par 2.3 *************************
    #INTEGER n,np,indx(n)
    #REAL a(np,np),b(n)
    #INTEGER i,ii,j,ll
    #REAL sum
    ii=0
    global amat,indx,rhs
    for i in range(1,n+1):
        ll=indx[i]
        sum=rhs[ll]
        rhs[ll]=rhs[i]
        if ii!=0:
            for j in range(ii,i):
                sum=sum-amat[i,j]*rhs[j]
        #hp     else if (sum.ne.0.) then
        elif abs(sum)>1.0e-30:
            ii=i
        rhs[i]=sum

    for i in range(n,0,-1):
        sum=rhs[i]
        for j in range(i+1,n+1):
            sum=sum-amat[i,j]*rhs[j]
        rhs[i]=sum/amat[i,i]
    #(C) Copr. 1986-92 Numerical Recipes Software Ea%].3.

def transp(c0, c, q, delta, nmax,ijmax, nc0, nci, ifrom, jto):
    # real c0(nmax), c(nmax), q(0:nmax, 0:nmax), delta(0:nmax, 0:nmax),
    # &     tran(ijmax)
    #  integer ifrom(ijmax), jto(ijmax)
    global tran

    for ij in range(1, nc0+1):
        j=jto[ij]
        tran[ij]=q[0,j]*(c0[j]+c[j])/2 + delta[0,j]*(c0[j]-c[j])
        tran[ij]=tran[ij]*360
    for ij in range(1+nc0,nc0+nci+1):
        i=ifrom[ij]
        j=jto[ij]
        tran[ij]=q[i,j]*(c[i]+c[j])/2 + delta[i,j]*(c[i]-c[j])
        tran[ij]=tran[ij]*360

def botd(tyear, nmax, n, ty, dt, s):
    # processing of the landsubsidance to the model results
    #
    global a,v
    for in_ in range(1,n+1):
        v[in_]=v[in_]-s[in_]*a[in_]*dt*lsub_v/(100*360)


##################
# INITIALISATION #
##################

# ########################################
# the following are static variables and #
# depending on individual tidal inlet    #
##########################################

nmax   = 3 # maximum array dimensions
rn     = 2 # power for equilibrium concentration
nodal_tide = 18.6 # years of nodal tide cycle
h      = 1.86  # tidal range !!!
anod   = 0.
anod   = 0.04*h # amplitude of nodal tide is 4%of mean tidal amplitude (dynamically changed in time loop)
ofile_vo = 'volume.out'   # output file volume change
ofile_ta = 'transport.out'   # output file tidal range
ofile_tr = 'tidal_range.out'   # output file tidal range

######################################################
# the following data is updated after file read ######
######################################################
ijmax=nmax
delta=numpy.zeros((nmax+1,nmax+1),  dtype=('f8')) #
dlta=numpy.zeros((ijmax+1),   dtype=('f8')) #horizontal exchange rate
a     = numpy.zeros((nmax+1), dtype=('f8')) # horizontal areas
v     = numpy.zeros((nmax), dtype=('f8')) # volumes
alfa  = numpy.zeros((nmax+1), dtype=('f8')) # equilibrium coeff.
ws    = numpy.zeros((nmax+1), dtype=('f8')) # vertical exchange velocity
trsp  = numpy.zeros((nmax+1), dtype=('f8')) # transport
tyear = []                                  # output years !!!
indx  = numpy.zeros((nmax+1),  dtype=('i4')) #
vcum  = numpy.zeros((nmax+1,10),  dtype=('f8')) # accumulated volume !!!
amat  = numpy.zeros((nmax+1,nmax+1),  dtype=('f8')) # needed in ludcmp and ...
ve    = numpy.zeros((nmax+1),dtype=('f8')) #
v     = numpy.zeros((nmax+1), dtype=('f8')) #
vv    = numpy.zeros((nmax+1),dtype=('f8')) # 
c0    = numpy.zeros((nmax+1),dtype=('f8')) #
tran  = numpy.zeros((ijmax+1),   dtype=('f8')) #
rhs   = numpy.zeros((nmax+1),     dtype=('f8')) #
c     = numpy.zeros((nmax+1), dtype=('f8')) #
ce    = numpy.zeros((nmax+1),dtype=('f8')) #
w     = numpy.zeros((nmax+1), dtype=('f8')) #
q     = numpy.zeros((nmax+1,nmax+1),   dtype=('f8')) #
en    = numpy.zeros((nmax+1),dtype=('f8')) #
ceglob= numpy.zeros((nmax+1),  dtype=('f8')) #
s     = numpy.zeros((nmax+1),     dtype=('i4')) #  
ifrom = numpy.zeros((ijmax+1),dtype=('i4')) # 
jto   = numpy.zeros((ijmax+1),  dtype=('i4')) #

h0     = 0.
p      = 0.    # in mkve calculated later on as vtidal
cce    = 0.0002   # overall or global equilibrium conc.
ystart = 2000. # year where sealevel increase starts
dt     = 36.   # timestep / days
year0  = 1970. # start year to establish equilibrium (updated by input file read later)
nt     = 5000. # number of time steps
yturn  = 2000. # start time of sea level rise
calls  = 1     # asmita time loop counter
trange = h     # tidal range
slr    = 0.    # later updated by input file
sed_incr_by_height_diff = 0 # increase equilibrium coefficient due to nodal tide higher mean sea level 
nodalt = 0

#############################
# read new input file #######
#############################
f=open('asmita.inp','r')
inp_data=[]
dummy=0

for l in f:
    l=l.rstrip() # remove empty spaces and newline
    if l[0]!='#':
       inp_data.append(l.split(','))
f.close()

# get asmita parameters
year0    = numpy.int(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '1'][0]][1])
yearl    = numpy.int(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '2'][0]][1]) 
dt       = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '3'][0]][1])
a[1]     = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '4'][0]][1])
a[2]     = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '5'][0]][1])
a[3]     = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '6'][0]][1])
v[1]     = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '7'][0]][1])
v[2]     = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '8'][0]][1])
v[3]     = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '9'][0]][1])
ws[1]    = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '10'][0]][1])
ws[2]    = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '11'][0]][1])
ws[3]    = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '12'][0]][1])
alfa[1]  = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '13'][0]][1])
alfa[2]  = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '14'][0]][1])
alfa[3]  = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '15'][0]][1])
cce      = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '16'][0]][1])
dlta[1]  = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '17'][0]][1])
dlta[2]  = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '18'][0]][1])
dlta[3]  = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '19'][0]][1])
nodalt   = numpy.int(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '20'][0]][1])
nt_amp   = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '21'][0]][1])
tr_inc   = numpy.int(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '22'][0]][1])
tr_inc_v = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '23'][0]][1])
slr_yes  = numpy.int(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '24'][0]][1])
slr_file = inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '25'][0]][1]
rn       = numpy.int(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '26'][0]][1])
h        = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '27'][0]][1])
lsub_yes = numpy.int(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '28'][0]][1])
lsub_v   = numpy.float(inp_data[[item for item in range(len(inp_data)) if inp_data[item][0] == '29'][0]][1])

print 'year0,yearl,dt,a[1],a[2],a[3],v[1],v[2],v[3],ws[1],ws[2],ws[3],alfa[1],alfa[2],alfa[3],cce,dlta[1],dlta[2],dlta[3],nodalt,nt_amp,tr_inc,tr_inc_v,slr_yes,slr_file,rn,h,lsub_yes,lsub_v',year0,yearl,dt,a[1],a[2],a[3],v[1],v[2],v[3],ws[1],ws[2],ws[3],alfa[1],alfa[2],alfa[3],cce,dlta[1],dlta[2],dlta[3],nodalt,nt_amp,tr_inc,tr_inc_v,slr_yes,slr_file,rn,h,lsub_yes,lsub_v

print 'check',v[1],v[2],v[3]
trange=h
nt=numpy.int(360*(yearl-year0)/dt)
print 'number of tme steps: ',nt
alfa_orig=[]
alfa_orig.append(-9999)
alfa_orig.append(alfa[1])
alfa_orig.append(alfa[2])
alfa_orig.append(alfa[3])

##################################
# read the landsubsidence data ###
##################################
# this is an array of years where data is available
# vcum is cumulative volume change due to landsubsidence
# pos means subsidence
#f=open(ifile2,'r')
#ctr=1
#tyear index starts with 1 hence:
#tyear.append(0)
#for l in f:
#    l=l.rstrip()
#    (t_year,vcum[1,ctr],vcum[2,ctr],vcum[3,ctr])=l.split(',')
#    tyear.append(float(t_year))
#    ctr+=1
#f.close()

##################################
# write initial volumes to file ##
##################################
f=open(ofile_vo,'w')
f.write(str(year0)+','+str(v[1]) +','+str(v[2])+','+str(v[3]))
f.close()

#####################################
# initialise trilateral project sea #
# level rise scenarios ##############
#####################################
#read file
if slr_yes:
    f=open(slr_file,'r')
    num_lines = sum(1 for line in open(slr_file))
    x=numpy.zeros((num_lines))
    y=numpy.zeros((num_lines))
    dummy=0
    for l in f:
        (x[dummy],y[dummy])=l.split(',')
        dummy+=1
    slr_scenario=interp1d(x, y, kind='cubic')
    f.close()

##################################################
############## MAIN LOOP #########################
##################################################

# set trold to actual tidal range of the tidal system
h0=h

#
# ASCII OUTPUT #
#
f_vo=open(ofile_vo,'w')
f_tr=open(ofile_tr,'w')
f_ta=open(ofile_ta,'w')

#
# NETCDF OUTPUT
#
nc_vol=netCDF4.Dataset('asmita.nc','w',format='NETCDF3_CLASSIC')
nc_vol.createDimension('time',None)
var_t=nc_vol.createVariable('time','d',('time',))
setattr(var_t,'units','seconds since ' + str(int(year0))  + '-01-01 00:00:00')
var_vol_flat=nc_vol.createVariable('vol_flat','d',('time',))
var_vol_channel=nc_vol.createVariable('vol_channel','d',('time',))
var_vol_delta=nc_vol.createVariable('vol_delta','d',('time',))
var_tra_flat=nc_vol.createVariable('transp_flat','d',('time',))
var_tra_channel=nc_vol.createVariable('transp_channel','d',('time',))
var_tra_delta=nc_vol.createVariable('transp_delta','d',('time',))
var_tr=nc_vol.createVariable('tidal_range','d',('time',))
var_slr=nc_vol.createVariable('mean_level_incr','d',('time',))

t_asm=[]
vol_flat=[]
vol_channel=[]
vol_delta=[]
tra_flat=[]
tra_channel=[]
tra_delta=[]
tr=[]
mwl=[]

first=0
v1_le_zero=1
v2_le_zero=1
v3_le_zero=1
for i in range(1,nt+1):
    if first and i%10==0:
        print t0,v[1],v[2],v[3]
        # fill NETCDF VARIABLE output
        t_asm.append(i*float(dt)*3600*24)
        vol_flat.append(v[1])
        vol_channel.append(v[2])
        vol_delta.append(v[3])
        tra_flat.append(trsp[1])
        tra_channel.append(trsp[2])
        tra_delta.append(trsp[3])
        tr.append(trange)
    first=1

    t0=(i-1)*float(dt)/360.+float(year0) # time value at n-1
    t1=t0+float(dt)/360.          # time value at n

    if nodalt: 
        anod=nt_amp*h/2. # amplitude of nodal tide is 4% of mean tidal amplitude 

        # maximum of nodal tide is ~1970, simulation starts at 1970, hence no 
        # phase shift 'fasey' needs to be aplied (fasey is in years)
        # vaue is read in input file fasey=8.

        # this is not 100% correct because the nodal tide maximum is not exactly 1970,
        # but until now...
        fasey=numpy.abs(years[0]-1970)
    else:
        anod=0.
        fasey=0.
    trange=h*(1.+float(anod)*numpy.cos(2*numpy.pi*(float(t0)-float(year0)-float(fasey))/nodal_tide))

        # increase equilibrium coefficients due to nodal tide increased sediment transport
        # by up to 20% and -20% 
        # A.P. Oost et al. (1993)
#        if sed_incr_by_height_diff: 
#            if trange-h>0:
#                percentage=0.7*(0.2*(trange-h)/anod)
#                alfa[1]=(1+percentage)*alfa_orig[1]
#                alfa[2]=(1+percentage)*alfa_orig[2]
#                alfa[3]=(1+percentage)*alfa_orig[3]
#            if trange-h<0:
#                percentage=0.7*(0.2*(trange-h)/anod)
#                alfa[1]=(1+percentage)*alfa_orig[1]
#                alfa[2]=(1+percentage)*alfa_orig[2]
#                alfa[3]=(1+percentage)*alfa_orig[3]

    if tr_inc: 
        # tidal range increase due to stronger increase of MHW than MLW, value per 
        # century given by variable trend in .inp file
        trange=float(trange)+float(tr_inc_v)*(float(t0)-float(year0))/100.

    # update slr by scenario 
    if slr_yes:
        if t0<2100:
            slr=slr_scenario(t1)-slr_scenario(t0)
        else: 
            slr=0.00001
    else:
        slr=0.

    if first and i%10==0:
        mwl.append(i*slr)

    # calculate transports, volumes

    asmita(dt,h,t0,t1,a,slr,cce,rn,p,tyear)

    h0=trange

    f_tr.write(str(t0)+','+str(trange)+'\n')
    f_ta.write(str(t1)+','+str(trsp[1])+','+str(trsp[2])+','+str(trsp[3])+'\n')
    f_vo.write(str(t1)+','+str(v[1])+','+str(v[2])+','+str(v[3])+'\n')

    if v[1]<0. and v1_le_zero:
        v1_le_zero=0
        print 'flat volume <0: ',t1,v[1]
        time.sleep(3)
    if v[2]<0. and v2_le_zero: 
        v2_le_zero=0
        print 'channel volume <0: ',t1,v[2]
        time.sleep(3)
    if v[3]<0. and v3_le_zero: 
        v3_le_zero=0
        print 'delta volume <0: ',t1,v[3]
        time.sleep(3)

    calls+=1


# finalise writing netcdf
print 'write to netcdf ...'
var_t[:]=t_asm[:]
var_vol_flat[:]=vol_flat[:]
var_vol_channel[:]=vol_channel[:]
var_vol_delta[:]=vol_delta[:]
var_tra_flat[:]=tra_flat[:]
var_tra_channel[:]=tra_channel[:]
var_tra_delta[:]=tra_delta[:]
var_tr[:]=tr[:]
var_slr[:]=mwl[:]
f_tr.close()
f_ta.close()
f_vo.close()
nc_vol.close()











